An Iterative Shrinkage Approach to 
Total- Variation Image Restoration 

Oleg Michailovich * 
October 26, 2009 

Abstract 

The problem of restoration of digital images from their degraded mea- 
surements plays a central role in a multitude of practically important 
applications. A particularly challenging instance of this problem oc- 
curs in the case when the degradation phenomenon is modeled by an 
ill-conditioned operator. In such a situation, the presence of noise makes 
it impossible to recover a valuable approximation of the image of inter- 
est without using some a priori information about its properties. Such a 
priori information - commonly referred to as simply priors - is essential 
for image restoration, rendering it stable and robust to noise. Moreover, 
using the priors makes the recovered images exhibit some plausible fea- 
tures of their original counterpart. Particularly, if the original image is 
known to be a piecewise smooth function, one of the standard priors used 
in this case is defined by the Rudin-Osher-Fatemi model, which results 
in total variation (TV) based image restoration. The current arsenal of 
algorithms for TV-based image restoration is vast. In the present paper, a 
different approach to the solution of the problem is proposed based on the 
method of iterative shrinkage (aka iterated thresholding) . In the proposed 
method, the TV-based image restoration is performed through a recur- 
sive application of two simple procedures, viz. linear filtering and soft 
thresholding. Therefore, the method can be identified as belonging to the 
group of first-order algorithms which are efficient in dealing with images 
of relatively large sizes. Another valuable feature of the proposed method 
consists in its working directly with the TV functional, rather then with 
its smoothed versions. Moreover, the method provides a single solution for 
both isotropic and anisotropic definitions of the TV functional, thereby 
establishing a useful connection between the two formulae. Finally, it is 
shown experimentally that, in the case when image degradation is caused 
by blur, the proposed method can provide restoration results of superior 
quality as compared to the case of sparse-wavelet deconvolution. 
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1 Introduction 



Both environmental effects and imperfections of image acquisition devices tend 
to degrade the quahty of imagery data, thereby making the problem of image 
restoration an integral part of modern imaging sciences. In particular, medi- 
cal imaging [1], astronomical [2] and laser [3] imaging, microscopy [4], remote 
sensing [5] , and photography [6] are all examples of applications in which the ne- 
cessity to improve the resolution and contrast of digital images routinely arises. 
Despite the relative "antiquity" of the theory of image restoration (with the 
first papers on the subject having been published as far back as at the end of 
the 60s [7,8]), there still exists a need for its further advancement via proposing 
new approaches as well as by improving the computational efficiency of existing 
ones. Addressing the second of the above two objectives forms the core of the 
developments presented here. 

The algorithm reported in this manuscript is based on a standard linear 
measurement model, in which the original image / and its measurements g are 
considered to be elements of a (finite dimensional) signal space U, and they axe 
assumed to be related according to 



where 7i : U — s- U is a bounded operator describing the effect of image degra- 
dation and e stands for both modeling and measurement noises. The problem 
of recovering / from g becomes particularly challenging in the case when the 
operator Ti. is either ill-conditioned or singular, in which case the problem of 
recovering / is commonly referred to as ill-posed [9]. On such conditions, a 
standard way to proceed with the solution of ([ij is to recover / approximately 
by solving the following variational problem 



where || • II2 denotes the ^2-norm and (/? : U — > M is a convex lower semicon- 
tinuous function on U whose role is to render the solution of Q unique and 
stable. Moreover, minimizing ip has the effect of restricting the solution of Q 
to functions of a predefined class which the original image / is believed to be- 
long to. In this case, the parameter A > (which is conventionally referred 
to as a regularization parameter) controls the balance between the model- and 
prior-dependent terms in (|2|. 

Probably the most renowned definition of ip as an ^2-type norm dates back 
to the works of A. N. Tikhonov [10]. In application to image processing, how- 
ever, this choice is rarely used in current practice because of the property of 
resulting solutions to have overly smoothed edges. In this regard, a more suc- 
cessful choice of ^p would be the one that allowed recovering / while maximally 
preserving its fine details. Following this line of considerations, Rudin et al [11] 
proposed to define the regularization functional ip to be the total variation (TV) 



g-H{/} + e, 



(1) 
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seminorm [12], viz. 



ip{u) = TY,{u) ||^|?^.P+|%|2||^, 



(3) 



where and Uy denote the partial derivatives of u G U and || • ||i stands for 
the £i-norm. The resulting model - commonly referred to as the Rudin-Osher- 
Fatemi (ROF) model - is nowadays considered to be one of the most fundamen- 
tal models of modern imaging sciences. Along with its isotropic version [13], 
i.e. 



the ROF model has proven to be an extremely useful tool in numerous ap- 
plications such as image de-blurring [2], "u -I- w" decomposition [14], super- 
resolution [15], and image impainting [16], just to name a few. 

Minimization of the TV-regularized functional E{f) = {5 ||^{/}^5ll2 + 
ATV(/)} (with TV e {TVj,TVa}) is known to be a relatively difficult opti- 
mization problem because of the non-differentiability of the TV regularizer. In 
order to overcome this difficulty, some methods substitute the absolute value 
function in ^ and Q by its smooth approximation [2, Section 8.2], [17,18]. 
Even though using such approximations provides an access to a variety of ef- 
ficient tools of smooth optimization, in order to perform stably, the resulting 
computational schemes require the use of proper preconditioning procedures, 
which makes these methods "costly" for fast processing of standard-size im- 
ages (e.g. 256 X 256 or 512 x 512). The same concern can be extended to the 
algorithms employing the tools of constrained optimization [19-21]. 

Practical problems related to the size of imagery data have motivated the 
community of imaging scientists to reconsider the potential of some first- order 
image reconstruction algorithms. In particular, the methods detailed in [22-25] 
are capable of finding a solution to ^ by means of simple recursive procedures, 
which make them particularly attractive for processing of large amounts of data. 
Unfortunately, these methods are only applicable to the de-noising setting (i.e. 
Ti. is an identity), with their extension to the case of non-trivial, rank-deficient 
H being currently considered impossible [25]. 

A different approach to the solution of ([2| was recently proposed in [26] 
based on the majorization-minimization (MM) method [27]. In this di- 
rect minimization of E{f) is substituted by recursively minimizing its quadratic 
majorizer whose minima can be found via solution of a system of linear equa- 
tions. For typical size images, however, the system can only be solved iteratively 
(using, e.g., the conjugate gradient algorithm), which substantially increases the 
overall computational cost of the procedure. It is also interesting to note that, 
even though derived from a different perspective, the method of [26] is essen- 
tially identical to the method of lagged-diffusivity [2,28]. Furthermore, a stable 
implementation of this method requires using a smoothed version of the TV 
functional (the fact not mentioned in [26]), which is necessary to prevent the 
diffusivity coefficients from becoming unbounded. 



ip{u) = TVa(u) := II |W: 




1' 



(4) 
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It is important to note that the recent interest in apphcation of MM-type 
strategies to image restoration seems to have been triggered by the works re- 
ported in [29,30] (see also [31] for a nice summary of the subject). In these 
works, the image restoration is performed under the assumption that / can be 
sparsely represented in the domain of a certain linear transform, which leads to 
the definition of 1^9 in ([2| as the £i-norm. The most remarkable result of these 
studies has been in showing that using the MM method allows one to solve the 
above problem by means of a simple first-order procedure. The latter - known 
as iterated shrinkage (aka iterative thresholding [32]) - consists of repetitive 
application of two simple steps: a back-projection correction and soft thresh- 
olding. Moreover, under a few standard assumptions, the iterated shrinkage is 
guaranteed to converge to a minimizer of the £1 -constrained cost functional. 

Despite the conceptual similarity between the TV- and £i-norm based reg- 
ularizers [33] , an iterative shrinkage approach to minimization of the cost func- 
tional in (|2| still seems to be missing. Thus, the main question addressed in 
the present study is whether or not it is possible to solve ^ for the case of 
non-trivial TC by means of an iterative shrinkage (IS) scheme. As will be shown 
below, in the discrete setting, the above question can be answered affirmatively. 
Accordingly, introducing an IS scheme for TV-based image restoration forms the 
main contribution of this work. Moreover, the proposed algorithm can be used 
to solve ([2]) for the cases of both tp — TV; and ip — TYa- In fact, the proposed 
numerical scheme will include a single scalar parameter which allows a transi- 
tion between the isotropic and anisotropic cases. Hence, another contribution 
of this paper consists in demonstrating a connection between TV^- and TVa- 
regularizers. Finally, it will be proven conceptually and experimentally that, 
in the case of ill-conditioned TC, the TV-based restoration can be expected to 
provide better reconstruction results as compared to the case of sparse wavelet 
regularization [30,31]. 

The remainder of the paper is organized as follows. Section II provides a 
number of essential technical details which are necessary for the developments 
in subsequent sections of the paper. An iterative shrinkage algorithm for the so- 
lution of Q with (p — TYi is detailed in Section III, whereas Section IV extends 
these results to the case of isotropic TV. Some important details regarding the 
implementation of the proposed method are discussed in Section V. Finally, the 
results of comparative experiments are summarized in Section VI, while Section 
VII finalizes the paper with a discussion and conclusions. 

2 Technical Preliminaries 

2.1 Signal Space 

In most of the practically important settings, images are finite dimensional 
objects. For this reason, we formulate our approach under the assumption that 
both original and measured images belong to the vector space of real-valued 
N X M matrices. Moreover, we endow this space (referred below to as U) with 
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the standard inner product (/, g) — ^n=o Sm=o^ fn,m gn,m and require that all 
the vectors in U are bounded and hence possess a finite f2-iiorm, defined in the 
standard way as ||/||2 = \/ (/, /)• Finally, the elements of the signal space are 
also constrained to have zero mean value, which leads to a formal definition of 
U as 

U={/eM^><*M (!,/>= 0, ||/||2<(X3}, (5) 

where 1 denotes an N x M matrix of ones. We note that the assumption of 
zero mean should not be regarded as a restrictive one, since in practice it is 
rarely a problem to subtract the mean value from a data image, as well as to 
re-normalize a zero-mean image to make its values saturate a required range, 
e.g., [0,255]. 

Let fx and fy be the partial differences of / taken in the column and row 
direction, respectively. Then, the discrete versions of the TV functionals (|3| 
and Q can be defined as 

TV.(/) = (l,^XmA^), (6) 

and 

TV„(/) = (1,|A| + |/,|). 

Consequently, the problem of TV-based image restoration of / 
restated as computing 

/tv = argmin , where 

fev 

E{f)=^-\\n{f}-g\\l + XTYif), (8) 
TVe {TV„TVa}, 

with the last expression in ([s]) meaning that the TV regularization can be either 
isotropic or anisotropic. 

In this paper we are particularly interested in the case when Ti. represents 
an operator of convolution, in which case the resulting restoration problem 
becomes that of image deconvolution. More specifically, H is assumed to be 
mean-preserving linear filtering, which implies = 1. Note that (subject 

to a proper normalization) the class of such blurs is relatively broad, including 
the important examples of moving-average, out-of-focus, and motion blurs, just 
to name a few. The mean-preserving property of H guarantees that z is 
the only vector contained in the intersection of the null spaces of ||7i{z}||2 : 
^NxM _^ ^ g^^^ TV(z) : M^xM M, which in turn suggests that, for any 
A > 0, the function ||H{z}||^ + ATV(z) : R^""^' ^ K is coercive and strictly 
convesj^ As a result, the optimization problem ^ is guaranteed to admit a 
unique global minimizer [25]. 

It should also be noted that, in the case oi g E V (which can always be en- 
forced by setting the mean value of the data image g to zero), the minimization 

^It is, in fact, a norm. 



(7) 

in ([!]) can be 
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in ([sj) can be performed over K^><^ (rather than U), with the solution guaran- 
teed to belong to U. This fact can be easily verified by contradiction as follows. 
Assume z e M^^^^ is a global minimizer of i? in ([s]). Denoting by /x 7^ the 
mean value of z (i.e. \i = {NM)~^ ^ ^n,m), the latter can be represented 
as z = z + ^ 1, where z G U. However, due to the orthogonality of U w.r.t. the 
subspace of constant images as well as due to the fact that TV(z) = TV(z) , it 
holds that 

E{z) = ^ \\n{~z} - + ^IIa^ III2 + TV(5) - E{z) + NMfi' > E{~z), 

which contradicts the assumption on z to be a global minimizer. 
The fact that, for g GU, the global minimizer 

/tv = argmin {E{f)} (9) 
belongs to U will play a key role in the derivations that follow. 



2.2 Gradient and Divergence Operators on U 

Let V denote the operator of discrete gradient defined in the standard manner 
as 

V : M^^*^ ^ (M^><*^)2 : / ^ ^/ = ( ) (10) 

for some standard definitions of the partial differences fx and fy (see below). 
The method proposed in this paper is based on the fact that the restriction of V 
to U is infective and hence invertible on its image. Let Vu denote the restriction 
of V to U and V := range(Vu) = range(V). Then, the above statement suggests 
that there exists a left inverse operator 

U:¥^V:v^(^ ^^u (11) 

such that 

U{Vf}^f, V/eU. (12) 

Below we explicitly construct the gradient V and the operator hi for two prac- 
tically important cases of replicative and periodic boundary conditions. 



2.2.1 Replicative boundary conditions 

In this case the discrete gradient V/ can be defined according to 

/yy /.\ J {fx)n,rn fn,rn /n— l,m; with /— l.rn /o,m f'io^ 

{vj)n,m ^ ] If \ _ r _ r -ih f — f ^ ^ 

\\Jy}n,rn Jn,m /n,m— 1; WILU Jn^~l /n,0 
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where n = 0, 1, . . . , — 1 and to = 0, 1, . . . , M — 1. Subsequently, congruent to 
the definition of (13 1, the discrete divergence operator div : (M^^^)^ ~* ^nxm 



can then be defined as 



div(v), where 

{Vx)l,ni, if 71 = 

(div(v))„ „ = <( (f:r)n+l,m " Mn.m, ifO<n<iV-l + (14) 
-{Vx)N~l,m, ifn = 7V-l 

{Vy)n,l; if TO = 

{Vy)n,rn+1 " (l'y)n,m, if < TO < A'/ - 1 
-{Vy)n,M-l' ifTO = Af-l 

for any v G (M^><*^)2. 

It is important to point out that the above definitions of the gradient and 
divergence operators are consistent with their continuous counterparts in the 
sense that —div constitutes the adjoint operator of V. Specifically, for any 
V e (M^><*^)2 and u e M^><*^ it holds that 

(V">v)(,j«xM)2 = (u,-div(v)) , (15) 

where {v,w) i^^nxm^ = {vx^w^) + {vy,Wy). 

Let VCT : M^^*^ j^AfxM ^y^q operator of 2-D discrete cosine transform 
(DCT) (as it can be implemented using, e.g., the dct2 function of MATLAB). 



Then, with the definitions ( |13[ ) and (14), it is straightforward to show that, for 
any u G R^^^ 

PCTjdiv (Vw)} = VCT{u} ■ W, (16) 

where the dot stands for element-wise matrix product and the elements of the 
N X AI matrix W are defined as 

P^,, = 2cos|^|+2cos|g|-4, (17) 

with fc = 0, 1,...,7V— 1 and ? = 0, 1, . . . , M — 1. It is important to point out 
that all the values Wk,i are strictly positive for fc + Z > 0, while Wo^ = 0. This 
fact makes it possible to define the integration filter Wi E M^^*^ as 

(PF,),,= |[2^°^{l^} + 2cos{f^}-4]-\ k + l>0 ^^g^ 



which satisfies 
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In combination with ( 16 ) and ( 19 1, the fact that, for any u G U, (X>CT{w})q q = 
suggests that 

u = VCT-^ {VCT {div(Vu)} • lyj , (20) 
which, in turn, leads to the definition of operator U as 

U : (M^^*^^ u . ^ ^ VCT-^ {I?Cr{div(v)} • W^} , 



(21) 



Note that 14 defined by (21 1 obviously satisfies (12 1. 



2.2.2 Periodic boundary conditions 



For the sake of completeness, we provide definitions analogous to (13), (14), and 
(20) for the case of periodic boundary conditions. Specifically, in this case, the 
gradient operator is defined as 

(V/) 



{fy)n,n 



fn,'. 
fn.'i 



with n = 0,1, . . . , N — 1 and m 
divergence operator is defined as 



fn~l,m, with /_i^„i — fN-l,m ^22^ 
/n,m-l, with fn,-l — fn,M~l 

- 0,1,...,M — 1, while the corresponding 



(div(v))„.„ 



{Vx)„+l^rn - {Vx)n,m, ifO<n<iV-l 
{yx)a,rn - {Vx)N~l,m, if U ^ N - 1 



iVy)r 



if < m < A/ - 1 



{Vv)n,0 - {Vy)n,M-l^ if m = M - 1 



(23) 
(24) 



for any v e (M^><^^)2. 

For the above definitions of the gradient and divergence operators it can be 



shown that, for any u G 



has 



VTT{div (Vm)} = VTT{u} ■ W, 



(25) 



where VTT : (^nxm jg ^-^^^ operator of 2-D discrete Fourier transform 

(DFT) (as it can be implemented using, e.g., the fft2 function of MATLAB) 
and 

(27:k\ 



W^fc.i =2cos. 



2 cos < 



(2t:1 



4, 



(26) 



A^- 1 and Z = 0,1,.. .,M- 1. 

the fact that {VTT{u})a^o 0,V'u G U 



2.2.1 



with fc = 0, 1, 

Similarly to the case of Section 
in conjunction with ( 25 1 leads us to conclude that 

u = VTT^^ {VTT {diy(yu)} ■ W,} 



with 



_ /[2cos{^} + 2cos{^}-4] \ k + l>0 



0, 



fc = / = 0. 



(27) 



(28) 
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holds for any image u in U. As a result, the operator U, defined as 

U : ^ U : V ^ VTT-^ {P^r{div(v)} • WJ , (29) 



satisfies the condition of (12 1 



2.3 Projection on the Range of V\ 



It is important to point out that the definitions oiU'm (21\ and ( 29 1 are nothing 
else but discrete counterparts of the well-known constructions used for solution 
of the Poisson equation in continuous variational analysis [34]. Furthermore, 
the discrete setting makes it straightforward to prove that, for an arbitrary 
V € (M^^*^)^, the vector u = Z^{v} constitutes a unique minimizer of the norm 
II Vu — v||(i{NxM)2 among all vectors of U, viz^ 



U{v} = argmin \ || Vm - v||L„xm)2 \ = 

uuhu, - v^Wl + \\uy - VyWl] , Vv e (]R^x^'^)2 (30) 



argmm • 

tiGU 



Consequently, ( 30 1 suggests that the composite operator \/lA defined as 

: (M^^*^)2 ^V: V{W{v}} (31) 



constitutes an orthogonal projection of (M^^*^)^ onto V. (The property (31 1 
of U, in fact, defines it as a left inverse.) It is also interesting to note that, 
considering v S {MJ^^^^)"^ to be a vector field, the operator VZ-/{v} sets to zero 
the rotational (solenoidal) component of v, and therefore VZ-/ is, in fact, an or- 
thogonal projector onto the subspace of irrotational (curl-free) vector fields [35]. 



3 Derivative Shrinkage 

Similarly to the case with many other image restoration problems, the solution 
of ([s]) can be given a statistical interpretation [36]. Particularly, from the view- 
point of Bayesian estimation, the isotropic regularizer TV; ([3| favors solutions / 
with independent and identically distributed values {y f)n,m- Moreover, at each 
pixel (n, m), the phase arctan(/j,//a;)„^m of {Vf)n,m is assumed to be uniformly 
distributed in (— tt, tt], while its magnitude |(V/)„_m| follows a Laplacian distri- 
bution. The anisotropic regularizer TVq (|4]), on the other hand, suggests that 
the partial differences fx and fx are mutually independent and identically dis- 
tributed according to a Laplacian law. Needless to say that the reconstructions 
corresponding to TV^ and TVq generally differ in their properties and appear- 
ance. In this paper, a unified solution will be given to address both above cases. 

■^Note that the uniqueness of u = V^{^f} as a minimizer follows from the equivalence of the 
norms ||V- ||(]gjvxAf)2 and || ■ ||2 in U due to the special structure of U (i.e. "no constant images 
are allowed" ) as well as because of U being a finite dimensional subspace. 
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For methodological reasons, however, the case of the anisotropic TV Q will be 
worked out first. 

Let f — {fx, fy)^ G V denote the gradient of /, and hence, with a slight 
abuse of nomenclature, one can say U{V f} = U{i}. With the use of this new 
notation, we first notice that 



|f||i-EK/-)' 



l(/.) 



TVa(/), 



(32) 



and, hence, the optimization problem ^ (with TV = TVa) can be replaced by 
an equivalent problem of the form 



rTv„ =argmm< -|lWMf}}-.g||2 + A||f||i 



(33) 



with the operator U defined by either (21 1 or (29). In other words, in (33 1 the 



minimization over / e U is replaced by minimization over its partial differences 
f G V. Note that the equivalence between the original problem and ( 33 1 is 



underpinned by the existence of lA that defines a one-to-one correspondence 
between V and the signal space U. It should be pointed out that, as long as the 



minimization domain in ( 33 ) is restricted to be V, the cost functional in ( 33 1 



remains coercive and strictly convex, which in turn guarantees the existence of 
a unique global minimizer in U. Moreover, frVa can be used to recover the 
solution /tVo to the original problem ([8]) according to 



/tv„ — ^{fxVa}- 



(34) 



The form of ( 33 1 can be additionally simplified via introducing the operator 



U as a composition of the convolution 7i and integration U opera- 



nt : \ 

tors, i.e. A{-} :— Ti{U{-}}. Using A allows (33 1 to be rewritten in a more 
standardized form as 



frVa = argmin {i?TV„ (f )} , where 
fev 

ET^A^) = \u{i}-9\\i + m\i- 



(35) 



The problem ( 35 1 has a format identical to that of the sparse-constrainted 



reconstruction problems of [29, 30] , and hence it can be solved by the method 
of iterative shrinkage. The fact, however, that Etv^ in (35 1 is minimized over 
V (rather than over {M.^^'^y) makes it necessary to supplement each step of 



the iterative shrinkage by the orthogonal projection onto V according to (31 1 



The resulting algorithm can be summarized as follows. Let Sr (with t > 0) be 
the operator of thresholding defined in the standard way as 



Sr{x) = sign(a;) (|x| - t)_ 



(36) 
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Then the proposed algorithm for iteratively solving ( 35 ) finds fTVa a station- 



ary point of the sequence of estimates produced by the following iterations 

f (*+^) - {f + -A* {.g - ^{f (*)}}} (37) 

f(*+i) = VZ^{f(*+5)}, 

where A* is the adjoint operator of A, and c is a positive scalar obeying c > 
The structure of (37 1 can be additionally simplified by introducing 
b := vA*{(7} and ??.(•) := In this case, the iteration procedure (37 1 

can be rewritten more concisely as 

f - VU {f + (b - 7^{f } } , (38) 



Note that (38), in fact, defines a map from V to itself, which can be shown 
to be non-expansive and asymptotically regular (see Remark 3.12 in [29]). In 
combination with the property of V being a convex and closed set, the above fact 



guarantees that the iterations in (38 1 converge to a minimizer of i?TVa I^Qi Prop. 
3.9]. 

To complete the description of the iterative shrinkage algorithm for TVa- 
based image restoration, the operators A* and TZ should be explicitly defined 
along with the constant c. To this end, we first recall that A is defined as a 
composition of the operators 7i and hi, viz. 

A:Y -.w [VTT-^ {PJ^T {div(v)} • Wi}] (39) 

The above operation can be implemented at the cost of an FFT-based convolu- 
tion if the operator 7i corresponds to periodic convolution. In this case, 7i can 
be represented by an x M matrix H of the DFT of its associated convolution 
kernel, which allows the operator A to be defined as 

A:N^V:^^ VTT-^ {DTT {div(v)} • A} , (40) 

with A being the frequency response of the composition of Ti and integration 
W^, i.e. 

A^W^-H. (41) 

It should be noted that periodic boundary conditions are common in image 
processing, since they allow substantially reducing the computational cost of 
filtering-type operations through the use of FFT. Moreover, there is a number 
of standard techniques which make it possible to adapt arbitrary images for the 
processing by periodic convolution [37]. For these reasons, the derivations below 
will be confined to the case of periodic convolution with the operator A defined 



by (401 and (41) 



To specify the adjoint operator A* of A, let v and u be two arbitrary vectors 
V and U, respectively. Then, 

{A{^},u) = {^,A*{u}) = (div(v),I?J-r-i {VTT{u} ■ A}) , (42) 
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and, therefore, 

A* -.V -.u^ -W {VTT-^ {DTT{u} ■ A]) , (43) 

where A denotes the complex conjugate of A. Note that, using the above defi- 
nition of A* , the vector b in ([38| can be (pre-) computed according to 



b = A*{g} = -V {DTT{g} ■ A]) , (44) 

in which case it appears to be unnecessary to preprocess g by setting to zero its 
mean value, as the "DC" component of'DJ-'T{g} is, in any event, multiplied by 
{A)ofl — 0. Moreover, by direct substitution one obtains 

7^ : V ^ V : V -V {VTT-^ {VTT {div(v) • }) , (45) 

with |Ap = A • A. Note that the computational cost of applying TZ is actually 
defined by the cost of one FFT-based convolution. 

Finally, to determine the range of admissible values of the parameter c in 



(38), we first note that, for an arbitrary u G U, it holds that 

-1 



A{A*{u}} = VTT-^ {VTT{u} ■ \A\^ ■ W) (46) 



with W given by ( 26 1 . This suggests that the composition AA* corresponds to 
convolution of an input image with VTl • W}. Consequently, using the 

fact that l^p -W = W,- \H\'^, one has 

\\AA*\\ -max(|W,| • |i?n„^„ <max(|T4^,|)„,„max(|i7n„_^, (47) 



Moreover, using the definition (28|, it is straightforward to show that 
max(|VKi|)„,m = 



2 - 2 cos 



27r 

max{A^, M] 



(48) 



Consequently, subject to the normalization max„.,„(|iJp)„^m — 1, it follows 
that the admissible values of c should obey c > [2 — 2cos (27r/ maxjA*", M})]~^. 

Algorithm 1 below provides an outline of the proposed method for TVa-based 
image restoration through iterative shrinkage. (Note that it is assumed that the 
blur operator has been normalized to have max„ ,„(|_ff p)„ .m — 1) The primary 
purpose of Algorithm 1 is to connect together the most important results on 
this section, while more general versions of the method will be discussed in the 
sections that follow. 



4 Extension to the Case of Isotropic TV 

4.1 Multidirectional Gradient 

The method of the previous section has been derived for the TV-regularizer in 
([8| equal to TVa as given by Q. It is known, however, that using the isotropic 
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Algorithm 1 TVa-based image restoration by iterative shrinkage 

c ^ [2 - 2 cos (27r/max{A^, M})]~^ + e (for some e > 0) 
h<=A*{g} (using (|44|) 
f ^ Vg (using ^) 
while "f keeps changing" do 

f 4= {f + c-i (b - 7^{f})| (using ^ and (Esl) 
f VU{{} (using ([22]) and ([29|) 
end while 

/tv„ ^ (using ([29|) 

Re-normahze /tVo (optional) 



TV-regularizer ^ can provide quahtatively different solutions to the restoration 
problem, ft is, therefore, tempting to extend the results of the preceding sections 
to the case of isotropic TV regularization, i.e. TV = TV^. 

To find a connection between TV^ and TV^, we take advantage of the fol- 
lowing identity 



it/2 

(|acos0 + bsm9\ + |6cos6' — asin 

9\)d9^Va^ + b^, (49) 



which holds for any pair of real numbers a and b. Using the fact that /J^^^ (cos 9 + sin 9) d9 = 
2, the equality (49 1 can be alternatively expressed as 

nTT/2 



J^' (|acos6» + 6sin0| + |6cos6l - asin6l|) d9 



J^^^ {cos 9 + sin 9) d9 



y/a^ + b^. (50) 



Now, let 9l — {9q, 9^, . . . , 9j^_^} be a set of L points uniformly distributed 
in [0,7r/2). In particular, we define 6*^ = 7rfc/2L, with /c = 0, 1, . . . , L - 1. These 
points can be used to construct a Riemannian approximation I{a,b;L) to (49 1 
(or, equivalently, to (50l) as given by 



J, ^. ^ Efc^o ( I « cos 9^ + b sin 9^\ + \b cos 9^ - a sin 9^ \ ) (7r/2L) ^ 

Et=oi^os9j^ + sin9^){7r/2L) 
EkZo {\acos9^ + b sin 9 j^\ + \b cos 9^ - a sin 9^ \ ) 

Due to the continuity of the integrands in (|50[), the Riemannian approximation 



is guaranteed to converge to Va^ + as L goes to infinity. Formally. 

lim /(a, b; L) = y/ a"^ + b^ . (52) 

L — )-oo 

Moreover, when considered as a function of (a, b) G M^, /(a, 6; L) represents an 
upper bound on \/ + for any value L > 1. This fact can be formalized in 
the following lemma. 
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Va2 + 62 I{a,b-L) 




Figure 1: The upper row of subplots compare the functions \/a^ + fe^ and 
I{a, b;L — 3) via visuaKzing them as surface plots. The lower row of subplots 
show the same functions as gray-scale images superimposed by their correspond- 
ing level set contours. 

Lemma 1 For any L > 1 and any (a, b) E M^; 

I{a, b; L) > ^/aP^P, (53) 

while the equality holds only for those (a, b) which satisfy either 

acosOj; + bsm9j^ ^0 or - asin6'^ + 6cos6'f' 0, (54) 

where 6^ = 7rfc/2L, with /c = 0, 1, . . . , L — 1. 

The proof of Lemma 1 is straightforward, yet technical. For the reason of 
space, therefore, the proof is omitted here. Instead, Fig. [l] compares the func- 
tions + and J(a, 6; L) (for the case L = 3) visualizing them as both 
surface plots (upper rows of subplots in Fig. [T]) and gray-scale images super- 
imposed by their corresponding level set contours (lower rows of subplots in 
the same figure.) One can see that, as suggested by the lemma, the functions 
y/a^ + b'^ and /(a, 6; L) coincide along the directions in M? which are defined by 
the angles {d^}\izl and {6*^ + 7r/2}^jrQ . It is also interesting to note that the 
level sets of /(a, h; L) provide a piecewise linear approximation to the level sets 
of \/a^~-\-W, which becomes progressively more accurate as L increases. 
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Using the approximating function /(a, b; L) allows an alternative definition 
of TV, as follows 



TV,(/) = (l, + /2) = ^lirn^ (l, L) 
which, in turn, leads to the following approximation of TV^ 



TV,(/)~^l,/(/„/^;L)^ 

L-1 

= (l , ^ I COS 0^ + /j, sin e^\ + \fy cos 61^ - sin 6^ 



(55) 



(56) 



fc=0 



with 



'L-1 

E 

.k=0 



-, -1 



COS! 



(57) 



The quality of the approximation ( 56 1 is supposed to become progressively bet- 
ter as L increases. In the experimental part of the paper, however, it will be 
shown experimentally that setting L — 3 results in image restoration practically 
indistinguishable from the case of exact TV^. It is also important to note that, 
for the case L ^ 1, one has 



L=l 



TV,(/). 



(58) 



In other words, setting L = 1 transforms the approximative TV functional into 
(56) into TVa- For the convenience of future referencing, let the approximative 
TV functional be denoted by TVl, viz. 

(59) 



TVi(/) = ^l,/(/„/,;L)^ 

L-1 

= (l , ^ I COS 0^ + /y sin e]:\ + \fy COS ei - sin 6^ 



fc=0 



Then, the principal properties of TVl can be summarized as 

TV,(/) < TVi(/) < TV„(/) 

Vi > 1, V/ e M^''*^ 



and 



TVl(/) 



TV„(/), ifL=l 
TV,(/), if L^(X3 



(60) 



(61) 



Note that (60 1 follows as a direct consequence of Lemma 1, while (61 1 is due to 
([52]). 

Finalizing this subsection, it is instructive to take a closer look at the struc- 
ture of the summands in the definition ( 59 1 of TVl , which have the form of 
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I /j; COS 0^ + /y sin 6*1' I + I cos — /a; sin 0^ I . Let V/, as previously, be the 
discrete gradient of some / e K^^*^ computed according to either ( |13| or (22 1. 



Now, let us rotate V/ through the angle dj^ in the counter-clockwise direction 

./ can be expressed via V/ and 9j^, and it is given 



Such a rotated gradient 

by 



fx COS ej^ 



fy cos - fx sin ( 



(62) 



Subsequently, given a set of L different angles 9l = {^fc}fc=o ' can compute 
the rotated gradients corresponding to each of the given 9^. Formally, we in- 
troduce the notion of multidirectional gradient (MDG) which is defined as the 
map given by 



t,NxM 



r,NxM\2L 



(63) 



Now, let f be the MDG of /, i.e. f = Vg^f = {v gL f,V gL f, . . . ,W ^lJ^ . 
Then, the ^i-norm of f in {E.^'^^'^Y^ can be defined in the standard manner as 



L-l 



iif II 1 = E E (I + «™ I + 1 A cos - s™ ^- 



(64) 



n,m k=0 



in which case ( 59 ) suggests that 

TVi(/) = dL||f||i. 



(65) 



The relation ( 65 1 establishes a connection between the TV^ functional and 



the ^i-norm which is necessary to derive a method for TV i-based image restora- 
tion by means of iterative shrinkage, which is detailed below. 



4.2 TV^-based Reconstruction via Iterative Shrinkage 

The fact that the MDG Vg^ depends linearly on V suggests that the null space 
of Vg^ consists of the subset of all constant images in M^^*^, and, therefore, 
both Vq^ and its restriction to U share the same range which we denote by 
Vl, i.e. Vl := range(Vg^). Note that by the definition of the range, for any 
V e Vl there exists u G U such that '^gj^u = v. Running a few steps forward, 
let us assume that the restriction of Vg^ to U is injective, and hence invertible 
on its image. In such a case, Wg^ has to have a left inverse ■ — > U whose 
defining properties are 

1. UL{VgJ} = f, for all /eU, 

2. For any v G {R^^'^'^f^, argmin„gu W^g^u - ^ UlM- 
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Consequently, using the operator Ul one can replace the original restoration 
problem ([s]) with TV — TVl by an equivalent problem 

fTv. - argmin ||H{Z^i{f}} - g\\l + [Xd^Wh] , (66) 

which appears to be in the format suitable for its solution via iterative shrink- 
age [29,30]. Note that given fxVi) the corresponding estimate /tVi, of / can 
be computed according to 

/TV, =Z^L{fTvJ, (67) 
which coincides with /xVa for L = \ and approaches /xVi when L increases. 



Therefore, by merely varying the value of L, (66 1 allows switching between the 
anisotropic and isotropic TV reconstructions. 

Throughout the rest of this section, we construct the operator IAl and pro- 
vide the definitions of operators analogous to the operators A, A* , and TZ of 
the preceding section. To proceed with the derivations, it is first necessary to 
find the divergence operator div^^ which is congruent with the definition of 
Vg . Such operator divg can be computed based on its property as an adjoint 
operator which requires that 

(V9^u,v)(RivxM)2i = (m, -divg^v) (68) 

holds for any u G M^^*^ and v e (E^^^-f)^^. Specifically, straightforward 
computations lead to the definition of divg^ as 

div,- : (R^-M^2L ^ jgATxM . ^ ^ ( gfe^o COS 9^ - v'^ sin 9^) \ 

(69) 

where div is given by either (14) or (23) in compliance with the correspond- 



ing definition of V, and the elements of v are assumed to be ordered as v 



(nfi „,o „.i „,i ii^^i v^^^^ 

yU^, U^, . . . , , Uy 



Given Vg^ and divg^, the operator Ul can be found based on its property as 
a projection operator. Specifically, let v be an arbitrary element in (M^^*^)^-'^. 

II 1 1 2 

Then, the vector w S U that minimizes the norm ||Vg^it — v|| ^jj„^jj,j^^2i should 
solve the systems of corresponding normal equations, which can be defined in 
the operator form as [38] 

div,-^ (Ve^u) =divgjv). (70) 
To solve the above system we will need the result of the following lemma. 



Lemma 2 Let L > I, and let and divg be the operators defined by (63) 



and (69) with respect to V and div given by either (13) and (14-) or (22) and 



(23), respectively. Then, for any w S U 

divg (Vg w) = idiv(Vw). (71) 
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Lemma 2 is proven via direct substitution of the definitions of V^^ and 
divg^ in the left-hand side of (71 1 with the use of some standard trigonometric 
equalities. Again, for the reason of space, the proof is omitted here. 

Using the result of Lemma 2 allows the normal equations ( 70 ) to be expressed 
in a different, yet equivalent form, viz. 

div(VM) = L"Miv0^(v). (72) 

Subsequently, depending on the type of boundary conditions, a u satisfying 
(72) can be found by means of either DCT or DFT transforms by virtue of 
the relations in (16) or (25). For the sake of concreteness, let the boundary 
conditions to be of the periodic type. In this case, the vector w G U uniquely 

II 1 1 2 

minimizing u - v , is found to be 



(73) 



where Wi is given by ( 28 1 . It should be emphasized that the uniqueness of the 
above solution is guaranteed by the fact that u is restricted to be an element of 
U. 



Based on (73), the left inverse operator Ul can now be defined as 



(I 



U : V 



1 



(74) 



Moreover, following the same line of considerations as in Section [23] the oper- 
ator of orthogonal projection from (M^^*^)^^ onto Yl — range(Vgj^) has the 
form of 



^S.i^L : (I 



i)JVxM\2L 



(75) 



As was already argued before, a considerable gain in computational efficiency 
becomes possible if the convolution operator H is defined to be periodic. In this 
case, one can define Al to be the composition operator given by 



A, 



V:^^n {UlM} = ^VJ^T-^ {VJ^T{divg^ (v)} • A} , (76) 



with A defined by (41). This allows the optimization problem (66) to be rede- 
fined in a more standardized way as 

1 



fTVi = argmin 



ALm-9\\2 + i^dL)\\{\\i 



(77) 



To find the solution fTVi of (77), the adjoint operator A}^ of Al and the 
composition TZl{-) ■— A*]^{Al{-}} need to be specified next. These operators 
can be shown to be respectively given by 



: U ^ Vi : M I 



--Vg^ {VTT{u} ■ A)) (78) 
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and 



7^L : ¥l ^ Vi : V 



(79) 



where, as previously, A stands for the complex conjugate of A and \A\'^ = A - A. 
It is important to note that, despite the use of L different gradient directions, 
the computation cost of applying TZ^ is still dominated by the cost of only one 
FFT-based convolution. 

Subsequently, given the above definitions and the theoretical guarantees 



of [29] , the solution to ( 77 1 can be found iteratively by the following recursion 

.{fW})}} (80) 



f - ^s^Ul [s^ {f 



(t) 



hr-n, 



where the soft thresholding St is defined by ( 36 ) , A*]^{g}, and c is required 

to satisfy c > It should be noted that the latter condition is crucial 

for the convergence of the iterative shrinkage procedure of (80 1. To determine a 
range of admissible values of c, we first notice that, for any arbitrary u G M.^^^ ^ 
one has 

{Al{u}} = -^VTT-^ {VTT{u} ■ W ■ \A\^} , (81) 
which, by comparison with (46 1, leads us to conclude that 

1 



WAlAIW = -\\AA* 



(82) 



Therefore, provided the convolution blur is normalized to obey max„ ,„(|_ff p)„^ 
1, and based on (48), one can conclude that c should be chosen to satisfy 

2tt 



1 



2 - 2 COS 



max{7V, M} 



(83) 



The central results of this section are summarized in Algorithm 2 below. 
Although the main purpose of Algorithm 2 it to establish connections between 
the main theoretical results of the paper, it can also be regarded as a "working 
prototype" of the proposed method. 



5 Technical Remarks 

5.1 MATLAB implementation 

For the sake of reproducibility of the results of this paper, some principal rou- 
tines needed for implementation of the proposed method are detailed next. In 
particular, in this subsection, we provide examples of MATLAB@ codes for 
computation of operators Vg^^, divg^, Ul, Al, A'^, TZl, and the update equa- 
tion ( [80| . (Note that the codes below have been optimized for clarity rather than 
for speed; substantial "speed-ups" are hence possible.) To this end, the vari- 
ables u and V will be used to denote generic elements of M^^*^ and (M^^*^)^-^, 
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Algorithm 2 TVi-based image restoration by iterative shrinkage 



Ol <= {7rfc/2L} 

Z^fe=0 



'L 

dL 
c ^ 

T ^ 

hL 

f 4: 



L-1 

k=0 



COS ( 



(1/L) [2- 2 cos (27r/max{A^,M})]" 

^Al{g} (using ([78])) 
Ve^g (using {fesl)) 



e (for some e > 0) 



while "f keeps changing" do 



Sr{i 



(bi 



f 

f ^ Vs^Ul{{} (using 
end while 

/tv, ^Z^L{f} (using ([74|) 
Re-normalize /tv^ (optional) 



7^L{f})} (using (361 and (791) 
631) and ([74])) 



respectively, with u handled as an N x M array, and v handled asanNxMx2*L 
array. 

Using the above notations, the MDG operator Wg^ can be computed by 
means of the m-function MDG which is given below. Note that the function 
computes the MDG using periodic boundary conditions. 

function [v] = MDG(u,L) 
theta = (pi/2/L)*(0:L-l) ; 
alpha=cos (theta) ; 
beta=sin(theta) ; 
[N,M] =size(u) ; 
ux=u-u([N,l:N-l] , :) ; 
uy=u-u(: , [M,1:M-1]) ; 
v=zeros(N,M,2*L) ; 
for k=l:L, 

v(: , : ,2*k-l)=alpha(k)*ux+beta(k)*uy; 

v(: , : ,2*k)=alpha(k)*uy-beta(k)*ux; 

end 

The multidirectional divergence (MDD) operator divg^ corresponding to Vg^ 
above can be implemented using the m-function MDD. 

function [u] = MDD(v,L) 
theta = (pi/2/L)*(0:L-l) ; 
alpha=cos (theta) ; 
beta=sin(theta) ; 
[N,M,K]=size(v) ; 
[ux,uy]=deal (zeros (N,M)) ; 
for k=l:L, 

ux=ux+(alpha(k)*v( : , : ,2*k-l)-beta(k)*v( : , : ,2*k)) ; 
uy=uy+(beta(k)*v( : , : ,2*k-l)+alpha(k)*v( : , : ,2*k)) ; 
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end 

u=(ux([2:N,l] , :)-ux)+(uy(: , [2:M,l])-uy) ; 



Having available the MDG and MDD functions, the computation of the remain- 
ing operators is straightforward. In particular, given the N x M matrices H and Wi 
of the transfer functions of the convolution and integration operators, respec- 
tively, and defining A = H. * Wi, the operator Al can now be computed using 
the m-function operator_A as given below. 

function [u] = operator_A(v, A,L) 
u=(l/L)*real(ifft2(fft2(MDD(v)) .*A))) ; 

It is important to note that the integration operator Ul can be implemented 
using the same m-function operator_A with the substitution of Wi for A. Finally, 
the m-function operator _A_star can be used to compute the adjoint operator 
Al 

function [v] = operator_A_star (u, A,L) 
v=(-l/L)*MDG(real(ifft2(fft2(u) .*conj (A)))) ; 

while the composite operator TZl can be implemented by the operator Jl func- 
tion defined as 

function [v] = operator_R(v, A,L) 
A2=A.*conj (A) ; 

v=(-l/L-2)*MDG(real(ifft2(fft2(MDD(v)) .*A2))) ; 

Using the above functions, the proposed method for TV-based image restora- 



tion can be implemented as a series of updates performed according to ( 80 1 . In 
particular, given the parameters c and tau defined by lines 3 and 4 of Algorithm 
2, respectively, and denoting by g and f the data image and its corresponding 
reconstruction, a total of K updates can be performed using the following simple 
codeEl 

b=operator_A_star (g,L) ; 

v=MDG(g,L); "/.initialization 
for k=l:K, 

v=wthresh(v+(l/c)*(b-operator_R(v,L)) , 's' ,tau) ; 
v=MDG(operator_A(v,Wi ,L) ,L) ; 

end 

f =operator_A(v,Wi,L) ; 

f =(255/range (f ( : ) ) ) * (f-min(f (:))); Zre-normalization 

Note that the last line of the above code forces the range of the recovered 
image f to fit the interval [0, 255], which by no means suggests it to be the only 
normalization scheme allowed. 



^Thc wthresh function is part of Wavelet ToolboxTM of MATLAB(1 
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5.2 Possible ways to improve the rate of convergence 



The computation of the multidirectional gradient and divergence operators can 
be performed with linear complexity as suggested by m-functions MDG and MDD 
above. Consequently, the computational cost of each update step of Algorithm 
2 is determined by the cost of 2-D FFT, and hence it has a complexity of 
O {NMlog{NM)). Although such a complexity can be considered as standard 
for many iterative methods of image deconvolution, it is important to spec- 
ify ways to reduce the total number of iterations performed by the algorithm, 
thereby minimizing its overall computational cost. To this end, it is first noted 
that the magnitude of the update term (1/c) (b^ — 7?,i,{f^*^}) and the threshold 
T = {XdLj/c in (80 1 are inversely proportional to the value of c. Therefore, the 
larger the value of c, the more significant is the change in f caused by each 
iteration. Thus, to maximize the effect of the update ( 80 1 , the value of c should 
be kept minimal. 

The requirement to minimize c appears to be at variance with the definition 



of its lower bound in (83 1, which suggests that c grows as max{A^, M} increases. 
This bound, however, should be considered as uniform in the sense of its being 
suitable for all iterations of the algorithm. It is known, on the other hand, that 
the method of iterative shrinkage belongs to the family of MM algorithms, in 
which a reduction in the value of an original cost functional is achieved through 
minimization of its local majorizers [27]. From this perspective, choosing c in 
accordance with ( 83 1 provides an a priori guarantee on all the majorizers to be 
convex. Alternatively, at a given iteration t, to get a reduction in the value of 



the cost functional in ( 66 1 , it is sufficient for c to satisfy 



Ct 



f _ f (t) 



> 



}:NxM\2L 



^^{f(t + l)_f(t)} 



(84) 



where the subscript t has been added to c to express its dependency on the 
specific iteration. Consequently, one can "fine-tune" the value of c based on 



( 84 ) by means of a simple back-tracking procedure as exemplified by Algorithm 



3. Note that, in this algorithm, the value of c is decreased exponentially (starting 
from an initial value obeying (83 1 ) by consecutively multiplying c by a reduction 
factor g (0, 1). 

An alternative way to minimize the value of c is to reduce the norm Hy^iyt^li 
through a sort of preconditioning. The latter can be defined by first noting that 
ll^i^J^II is determined by the maximum value of Wi ■ \H\'^ . Moreover, since 
in most of the practically important cases, the frequency response H can be 
normalized to satisfy < {\H\'^)n,m < 1, it is mainly the values of Wi which 
dominate the maximum of Wi ■ \H\'^ . Being an integration filter in nature, Wi 
tends to amplify the lower frequencies. As a matter of fact, the values of Wi 
are sharply peaked in a neighborhood of the zero frequency, as shown by Fig. [2] 
The figure visualizes Wi as a gray-scale image for the case oi N = M = 256. 
One can see that the overwhelming portion of the values of Wi are relatively 
small (as indicated by the black color), being equal to approximately 0.125. On 
the other hand, the values of the white pixels (which are few and hardly visible 
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Algorithm 3 Local adjustment of the value of c via back-tracking 
1: c ^ (1/L) [2 - 2cos(27r/max{A^, Af})]"^ + e (for some e > 0) 



T 4= {\dL)lc 

temp ^Sr{i+ {hL - c-l7^L{f }) } 
r <^ temp — f 

while c||r||^jj„xM)2L > (7^L{r}, r)(jjjvxAf )2l 
c lie (for some /i e (0, 1)) 
r ^ {XdL)/c 



do 



temp 5^ |f + (b/^ 

r <^ temp — f 
end while 
f Ve^Z^L{temp} 



7^L{f})} 



Integration filter Wi 




Figure 2: The frequency response of the integration filter Wi for the case of 
N = M = 256. 



without a proper zoom) are close to 1160.1, which causes the norm Hy^Ly^^li 
to have a relatively large value of 553.38, as computed according to (48 1 and 



(82) for L = 3. Consequently, the large values of are, in fact, due to a 



negligibly small number of high-amplitude values of Wi located around the DC. 

The property of |if p of being "flat" around the zero frequency implies that 
its values cannot "counterbalance" the high amplitudes of Wi. On the other 
hand, one can pre-convolve the data image g with another auxiliary filter whose 
DFT He is designed so that the product Wi ■ |-ffcP has a smaller maximum 
value than that of Wi. In this case, the resulting operator norm will 
be defined by the maximum value of Wi ■ |i?cP • \H\'^, which can be set to be 
reasonably small. It goes without saying that such a preconditioning changes 
neither the format of the restoration problem nor of its solution except for the 
obvious need to replace the original "blur" H hy a new one, i.e. H ■ He. 



Additional techniques for further speeding-up the convergence of (80) can 
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be borrowed from the field of projection methods of convex optimization [39]. 
In particular, let f ^*-* and f be two successive outcomes of the shrinkage 
procedure (80 1. Then, the next estimate of f can be found as a minimizer 
of the TVi-functional along the direction d(*+i) = f(*+i) - f(*). Note that, 
provided both f and f(*+i) are in V^, all solutions of the form f + ad(*+^) 
(with a S K) belong to Yl as well, and therefore the above line search can be 
performed without using intermediate projection steps. It was demonstrated 
in [40], that augmenting the shrinkage operation by the line search can result 
in substantial increase in the rate of convergence of iterative shrinkage. 

Finally, we note that the iterative shrinkage ( 80 1 belongs to the family of so- 
called one-step shrinkage schemes. It was recently shown in [41], that the latter 
can be extended to iwo-stage shrinkage schemes, which have considerably higher 
rate of convergence and, therefore, provide additional means to further speed 
up the implementation of the proposed method for TV-based image restoration. 



6 Results 

6.1 Reference methods 

In this section, the theoretical results derived in the preceding sections of this 
paper are supported by a number of experimental results. In particular, we 
first show that image restoration by means of the proposed TV-based iterative 
shrinkage (TVIS) method with TV = TVl and L — 3 provides restoration 
results virtually indistinguishable from the results obtained using alternative 
methods of solving ([s]) with TV = TV^. Additionally, it will be shown that, 
in the case of strong convolutional blurs, the TVIS method can provide more 
valuable restoration results as compared with the method of [30], which will 
be referred below to as the sparse wavelet iterative shrinkage (SWIS) method. 
Some minimal details about the references methods are given next for the sake 
of presentational completeness. 

6.1.1 The method of lagged diffusivity 

This method (which seems to have been first proposed by C. Vogel in [28]) 
is known to be one of the standard approaches to the solution of ([s]) with 
TV = TVi. The method is based on the first-order optimality condition for 
E{f) in ([s]), which is given by 

H*{^{/}-5}-Adiv(^^) =0, (85) 



where ||V/|| = y/^ + fy - Consequently, the global minimizer of E{f) is found 
as a stationary point of a sequence of solutions to the following system of equa- 
tions 



{w{/(*+^)}}-Adiv(^^^)=H*{.9} (86) 
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solved w.r.t. Z*-*"''"'^-'. Note that in (86 1, the result of previous iteration J*^*) is 
considered to be constant ("frozen"), which makes (86) be a linear operator 
equation that can be solved iteratively by means of, e.g., conjugate gradient 
algorithm. It is also worthwhile noting that the method of [28] was recently 
rediscovered in [26], where the same iterative procedure (86 1 is derived using 
the majorization-minimization (MM) technique. What appears to be omitted 
in [26], however, is mentioning the fact that a stable implementation of (86 1 
requires replacing the absolute value ||V/'*''|| by its strictly positive approxima- 



tion ||V/(*)|| w V(/i*^^^)2 + (/i*+^V + e, for some < e < 1. Thus, strictly 
speaking, as long as e > 0, the global minimizer of E{f) in ([s]) and the station- 
ary point of (86 1 cannot be guaranteed to be identical, in general. It is possible, 
however, to converge to a close vicinity of the true global minimizer using a 
smooth relaxation procedure in which a stationary point of (86 1 is found for 
some value of e > 0, followed by decreasing e (e.g., e <J= e/2) and, subsequently, 
computing a new stationary point of ( 86 1 for the new e, while using the previous 
stationary point as an initialization. In the present paper, the above relaxation 
"cycles" were performed to reduce the value of e from 10~^ to 10~^. 

Finally, we note that both the method of lagged diffusivity [28] and the 
proposed TVIS method require a definition of the regularization parameter A. 
An optimal value of A can be elicited based on the theory of Bayesian estimation, 
according to which the TV-based image restoration assumes (||V/||)„^ to be 
i.i.d. Laplacian, namely 

(l|V/||)„,,„ ^exp| 1, (87) 

where /3 > is a scale parameter of the distribution which can be estimated as 
(3 « y/o.S crji^ , with fyv/ii being the sample variance of |1V/||. Therefore, if 

the variance of the additive noise in ([ij is equal to cr^, then from the viewpoint 
of MAP estimation, the optimal A should be set to be equal to A = cr^//?, and 
this is how it was done in the present study. 



6.1.2 Sparse wavelet iterative shrinkage (SWIS) 

Let {V'fcjfcer be a dense set of wavelet functions in M^^*^, and ^' be their 
associated synthesis operator defined as 

feer 

In the present study, vp is defined to correspond to a stationary separable wavelet 
transform, which implies that the set {ipk}ker is overcomplete, and hence there 
is no unique way to represent an arbitrary / S M.^^^ in terms of ipk- However, 
if / can be sparsely represented by such a wavelet dictionary, then a useful 
approximation of / in ([T|) can be computed as 

argmin {l\\H{'f{x}}~g\\l + X\\xh], (89) 
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where A > is a regularization parameter analogous to that in ([8|. Moreover, 
it was proven in [29,30] that the above minimization problem can be solved via 
iterative shrinkage performed according to 

x^'^'^ = S,/c + {g - , (90) 

where .4{-} is the composition of H and "if, A* is the adjoint of A, Sr is given 
by (36 1, and the constant c satisfies c > 

Note that, from the perspective of Bayesian estimation, the representation 
coeSicients x in ( [89| are assumed to be i.i.d. Laplacian, which justifies setting 
A to be equal to ay /3, with /3 being the scale parameter of the distribution of x. 
To estimate f3, all the test images used in this study were first processed by the 
Basis Pursuit algorithm [42] to find their corresponding sparse representations. 
Subsequently, the sample variances of the sparse coeSicients were estimated, 
followed by estimating the respective f3 as y^O.ba'^. 

6.2 TVIS versus the method of lagged diffusivity 

Practical implementation of the method of lagged diffusivity [28] requires one 
to optimally preset a number of "internal" algorithmic parameters such as a 
convergence rule for the linear solver, the rate of relaxation for e, etc. Besides the 
somewhat esoteric nature of this preset, the latter also appears to be problem- 
dependent, which implies the necessity to adjust the algorithmic parameters 
for different data sets. The proposed TVIS algorithm, on the other hand, has a 
very simple structure which is devoid of any parameters that could be potentially 



dependent on data. Moreover, in Section 5.2 a number of means to increase the 
rate of convergence of the TVIS algorithm were proposed, which could be useful, 
for example, in the case of processing large sets of imagery data. We reserve a 
thorough discussion of that and related matters for another occasion. Instead, 
in this section, we focus on finding the minimal number of multidirectional 



gradients L in ( 56 1 for which the restoration results obtained with the reference 
method of [28] and TVIS can be considered as comparable. 

In the examples below, two different values of the parameter L of the TVIS 
algorithm are used, namely L ^ I and L = 3. Note that, in the case when 
L = 1, TVIS performs anisotropic TV-based image restoration, whereas it ap- 
proximates the solution of the isotropic TV-based image restoration problem 
when L = 3. (The above two settings will be referred below to as TVIS-1 and 
TVIS-3, respectively). The legitimacy and usefulness of the above approxima- 
tion are demonstrated through our first example in Fig. [3] The upper row of 
subplots in the figure shows an image of glomerulus, its blurred version, as well 
as the noisy image obtained by contaminating the blurred image with a white 
Gaussian noise. The blurring artifact was modeled by convolving the test image 
with a Gaussian kernel of standard deviation 0.8. It should be noted that, in this 
case, the condition number of the corresponding convolution operator Ti. (which 
can be computed as the ratio of the maximum and minimum values of the DFT 
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Zoom A Zoom B Zoom C 




Figure 3: (Upper row of subplots) Original image of glomerulus, its blurred and 
noisy versions (PSNR — 21.3 dB); (Middle row of subplots) Results of image 
restoration by (from left to right) TVIS-1, TVIS-3, and MLD; (Lower row of 
subplots) Zoomed segments of the estimated images as indicated by the yellow 
rectangles. 
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TVIS-1 estimate TVIS-3 estimate MLD estimate 




Figure 4: (Upper row of subplots) Original image of peppers, its blurred and 
noisy versions (PSNR = 24.5 dB); (Middle row of subplots) Results of image 
restoration by (from left to right) TVIS-1, TVIS-3, and MLD; (Lower row of 
subplots) Zoomed segments of the estimated images as indicated by the yellow 
rectangles. 
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of the Gaussian kernel) is equal to 138.4, which is relatively small. The blur 
artifact in Fig. [3] therefore, can be classified as mild. The noise contamination, 
on the other hand, is relatively strong, resulting in the peak signal-to-noise ratio 
(PSNR) of 21.3 dB. 

The results of TV-based restoration of the image of glomerulus by TVIS-1, 
TVIS-3, and the method of lagged diffusivity (MLD) are shown in the mid- 
dle row of subplots in Fig. [3] For all these methods, the same regularization 
parameter A = 0.028 was used, which was found based on the methodology de- 
tailed in Section |6.1| Interestingly enough, all the estimated images appear to 
be virtually indistinguishable. A closer examination, however, reveals that the 
anisotropic restoration by TVIS-1 has more "squarish" details as compared to 
the restorations by TVIS-3 and MLD. (An example of such a behavior is demon- 
strated by the zoomed segments of the image estimates shown in the lower row 
of subplots in Fig. [3j) The restorations obtained with TVIS-3 and MLD, on 
the other hand, are very close to each other, with a relative error between them 
being less than 0.5%. Moreover, both methods resulted in the same PSNR of 
24.4 dB. 

Despite the profound similarity between the restoration results in Fig. [3j 
the conclusion that TVa- and TV^-based image restorations are similar would 
be rather premature as proven by the example in Fig. |4] (whose composition is 
analogous to that of Fig. [3|. In this case, the standard image of peppers has 
been blurred by a Gaussian kernel of standard deviation 1.2, with the condition 
number of its associated H being equal to 326217.8. The addition of white 
Gaussian noises gave rise to the PSNR of 24.5 dB. 

Analyzing the restoration results shown in the middle row of subplots in 
Fig. |4] one can hardly see any difference between the estimates provided by 
TVIS-3 and MLD. Both methods resulted in the same improvement of 4.3 dB 
in terms of the PSNR, which is natural considering the fact the relative error 
between the TVIS-3 and MLD restorations was found to be below 0.5%. The 
anisotropic nature of TVIS-1, on the other hand, can now be clearly seen. This 
nature manifests itself in the irregular behavior of the image contours, which 
are smooth in the original scene. In this example, therefore, the restoration by 
TVIS-3 should be preferred over that by TVIS-1. 

6.3 TVIS versus SWIS 

Based on the results of the preceding subsection one can conclude that, for 
the case of L = 3, the estimations obtained with the TVIS algorithm are very 
close to what can be obtained using alternative methods for TV^-based image 
restoration. In this subsection, therefore, only the performance of TVIS-3 is 
compared with that of the method of image deconvolution by sparse wavelet 
iterative shrinkage (SWIS) [30,31]. To sparsely represent test images, stationary 
(aka non-decimated) separable wavelet transforms with three resolution levels 
were used. Note that, in this case, the number of wavelet coefficients is ten 
times higher than the number of image pixels, which suggests that the wavelet 
transform has the overcomplete factor of 10:1. The TVIS-3 algorithm, on the 
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Figure 5: (Subplot A) Histogram of the stationary wavelet coefficients of "Cam- 
eraman" computed by BP with the "Symlet-1" wavelet; (Subplot B) Histogram 
of the stationary wavelet coefficients of "Cameraman" computed by BP with the 
"Symlet-4" wavelet; (Subplot C) Histogram of the stationary wavelet coefficients 
of "Shepp-Logan" computed by BP with the "Symlet-1" wavelet; (Subplot D) 
Histogram of the stationary wavelet coefficients of "Shepp-Logan" computed by 
BP with the "Symlet-4" wavelet. 



other hand, has smaller storage requirements, as it operates with a total of six 
partial derivatives, resulting in the overcomplete factor of 6:1. 

At the heart of the SWIS method is the assumption that the image of in- 
terest can be sparsely represented in the domain of a wavelet transform. To 
verify this assumption, the basis pursuit (BP) algorithm [42] was employed to 
find sparse representations of the standard "Cameraman" and "Shepp-Logan 
phantom" images in the domain of the wavelet transforms corresponding to 
the nearly-symmetric wavelet of I. Daubechies having 1 and 4 vanishing mo- 
ment^ (Since in the standard nomenclature, the above wavelet functions are 
referred to as Symletl and Symlet4, the corresponding SWIS algorithms will 
be referred below as SWIS-SYMl and SWIS-SYM4, respectively). Subplot A 
and B of Fig. |5]show the histograms of the wavelet coefficients of "Cameraman" 
corresponding to Symletl and Symlet4, respectively, while the histograms of 
the wavelet coefficients of "Shepp-Logan phantom" computed using the above 
wavelets are shown in Subplots C and D of the same figure. The profoundly 
super-Gaussain appearances of all the histograms suggests that both test images 
are indeed sparsely represent able by the wavelet transforms in use. 

The upper row of subplots of Fig.[6]show the original image of "Cameraman" 
as well as its blurred and noisy versions. In this example, the "mild" Gaussian 

*In the case of one vanishing moment, the wavelet is identical to the "Haar" wavelet. 
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Figure 6: (Upper row of subplots) Original image of "Cameraman", its blurred 
and noisy versions (PSNR = 22.4 dB); (Middle row of subplots) Results of im- 
age restoration by (from left to right) SWIS-SYMl, SWIS-SYM4, and TVIS-3; 
(Lower row of subplots) Zoomed segments of the estimated images as indicated 
by the yellow rectangles. 
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Figure 7: (Upper row of subplots) Original image of " Shepp-Logan phantom" , 
its blurred and noisy versions (PSNR = 19.0 dB); (Middle row of subplots) 
Results of image restoration by (from left to right) SWIS-SYMl, SWIS-SYM4, 
and TVIS-3; (Lower row of subplots) Zoomed segments of the estimated images 
as indicated by the yellow rectangles. 
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blur of standard deviation 0.8 was used to smooth the image, followed by addi- 
tion of while Gaussian noise giving rise to the PSNR of 22.4 dB. Subsequently, 
the image was subjected to the restoration procedures of SWIS-SYMl, SWIS- 
SYM4, and TVIS-3 algorithms, whose results are depicted in the middle row of 
subplots in Fig. [6] In all the above cases, the regularization parameter A was 
computed according to the methodology discussed in Section |6.1| Specifically, 
for SWIS-SYMl, SWIS-SYM4, and TVIS-3, the optimal values of A were found 
to be equal to 0.1335, 0.1505, and 0.024, respectively. 

Analyzing Fig. [6] one can see that the best restoration results in terms of con- 
trast improvement and noise reduction are obtained by TVIS-3, with its PSNR 
equal to 26.8 dB. At the same time, both SWIS-SYMl (PSNR = 22.8 dB) and 
SWIS-SYM4 (PSNR = 22.4 dB) suppress the additive noise, while being unable 
to surmount the effect of blur. It is also interesting to note that SWIS-SYMl 
provides a "sharper" reconstruction as compared with SWIS-SYM4, while the 
latter results in smoother estimation of the uniform background. 

The inability of SWIS to overcome the effect of strong blurs is not occa- 
sional. This fact is further supported by Fig. [7] which demonstrates the results 
of restoration for "Shepp-Logan phantom" . In this example, the (relatively 
strong) Gaussian blur of standard deviation 1.2 was used to smooth the orig- 
inal image, which together with noise contamination resulted in the PSNR of 
19.0 dB. Analyzing Fig. (Tj one can once again see that TVIS-3 is capable of 
accurately recovering the piecewise constant structure of "Shepp-Logan phan- 
tom" (PSNR = 23.9 dB), while neither SWIS-SYMl (PSNR = 20.8 dB) nor 
SWIS-SYM4 (PSNR = 20.9 dB) could de-noise and sharpen the image to a 
similar extent. Moreover, we see again that SWIS-SYMl produces "sharper" 
(yet noisier) restoration as compared to SWIS-SYM4. 

The reason why SWIS cannot reduce strong blurring artifacts does not lie 
in the nature of the particular solution method used, but is intrinsic in the 
way the regularization is performed. Specifically, for the sake of illustration, 
let us assume that 7i in ([T]) is a low-pass half-band filter and no noise has been 
observed in the measurement of 7i{/}, i.e. e = 0. In this case, there is an 
infinite number of possible solutions to 7i{/} — g, all of which are different 
over the null space of H containing high-frequency signals whose spectra vanish 
below the cut-off of tt/2. Among these candidate solutions, we are interested 
in specific two, namely / and g (note that Ti.{g} = g, since Ti. is nilpotent). 
The sparse regularization via ^i-norm minimization can "prefer" / over g only 
if II /111 < \\g\\i- This is, however, not true in general. For example, in the case 
of Shannon multiresolution [43, Ch.VII], all the wavelet coefficients of / and g 
would be identical, except for the highest resolution level, where the coefficients 
of g would necessarily vanish as opposed to those of /. Obviously, in this case, 
ll/lli ^ \\g\\i holds for any pair {f,g), and, therefore, the minimization of the 
^i-norm would result in 5 as a final solution, not /. 

Even though ideal blurs are rare in practice (as well as the use of Shannon 
wavelets), the above considerations remain relevant in the case of relatively 
strong blurs and regular wavelet analysis. This seems to be the reason for 
which the examples in [30] and [31] used relatively weak blurs and irregular 
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Haar wavelets. In the case of heavier blurring artifacts, however, the assumption 
on the sparsity of wavelet coefficients (as measured by the ^i-norm) does not 
seem to be sufficiently strong to guarantee a successful image restoration as 
exemplified by Fig. [6] and Fig. [7] Under such conditions, performing TV-based 
image restoration (by means of, e.g., TVIS) should be considered as a more 
effective alternative. 



7 Discussion and Conclusions 

Among the existing methodologies for restoration of digital images, the TV- 
based method of [11] is known to be one of the most fundamental techniques. 
In this method, the instability of image restoration caused by the property of 
7i in {[ij being a poorly conditioned operator is overcome by using a priori 
information. The latter is introduced in the form of a requirement on the 
recovered image to have a small total-variation norm. In general, the solution 
of such a regularized problem depends on the definition of the TV-norm, which 
can be either anisotropic or isotropic. Thus, for some images, both definitions 
can result in very close estimations, while for others, the results may differ 
dramatically. In this respect, the proposed TVIS method is advantageous, for 
it allows one to switch between the above two settings by merely changing 



the value of L in (56 1. Moreover, it was demonstrated both conceptually and 
experimentally that, for L — the TVIS algorithm provides estimation results 
which are virtually indistinguishable from the results obtained by alternative 
methods of TVi-based image restoration. 

Yet another useful feature of the TVIS algorithm consists in its particularly 
simple structure, which requites only a recursive application of linear filtering 
and soft-thresholding. Furthermore, as opposed to many alternative methods of 
TV-based image restoration, TVIS does not require presetting any algorithmic 
parameters that could potentially depend on the properties of imagery data. 

The rate of convergence of the TVIS algorithm is defined by the image size as 



well as by the value of L according to ( 83 ) . In Section 5.2 however, a number of 



practical ways were detailed, using which one can substantially reduce the total 
number of iterations required by the algorithm. For the reason of limited space, 
a comparative analysis of these "speedup" schemes has not been included in the 
present paper, and it will be published elsewhere. However, even in its current 
form (see Algorithm 2), the proposed method has a complexity compared to 
that of the sparse wavelet iterative shrinkage in [29-31]. Yet, as opposed to 
SWTS, TVIS is capable of rehably recovering the images corrupted by relatively 
strong noises and blurring artifacts, as demonstrated by Fig. |6]and Fig. [7] 
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